<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">


<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
    
    <title>Chapter 7 BLAST &mdash; Biopython_en 1.0 documentation</title>
    
    <link rel="stylesheet" href="_static/default.css" type="text/css" />
    <link rel="stylesheet" href="_static/pygments.css" type="text/css" />
    
    <script type="text/javascript">
      var DOCUMENTATION_OPTIONS = {
        URL_ROOT:    './',
        VERSION:     '1.0',
        COLLAPSE_INDEX: false,
        FILE_SUFFIX: '.html',
        HAS_SOURCE:  true
      };
    </script>
    <script type="text/javascript" src="_static/jquery.js"></script>
    <script type="text/javascript" src="_static/underscore.js"></script>
    <script type="text/javascript" src="_static/doctools.js"></script>
    <link rel="top" title="Biopython_en 1.0 documentation" href="index.html" />
    <link rel="next" title="第8章 BLAST和其他序列搜索工具(实验性质的代码)" href="chr8.html" />
    <link rel="prev" title="Chapter 6 Multiple Sequence Alignment objects" href="chr6.html" /> 
  </head>
  <body>
    <div class="related">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="genindex.html" title="General Index"
             accesskey="I">index</a></li>
        <li class="right" >
          <a href="chr8.html" title="第8章 BLAST和其他序列搜索工具(实验性质的代码)"
             accesskey="N">next</a> |</li>
        <li class="right" >
          <a href="chr6.html" title="Chapter 6 Multiple Sequence Alignment objects"
             accesskey="P">previous</a> |</li>
        <li><a href="index.html">Biopython_en 1.0 documentation</a> &raquo;</li> 
      </ul>
    </div>  

    <div class="document">
      <div class="documentwrapper">
        <div class="bodywrapper">
          <div class="body">
            
  <div class="section" id="chapter-7-blast">
<h1>Chapter 7  BLAST<a class="headerlink" href="#chapter-7-blast" title="Permalink to this headline">¶</a></h1>
<p>Hey, everybody loves BLAST right? I mean, geez, how can get it get any
easier to do comparisons between one of your sequences and every other
sequence in the known world? But, of course, this section isn’t about
how cool BLAST is, since we already know that. It is about the problem
with BLAST – it can be really difficult to deal with the volume of data
generated by large runs, and to automate BLAST runs in general.</p>
<p>Fortunately, the Biopython folks know this only too well, so they’ve
developed lots of tools for dealing with BLAST and making things much
easier. This section details how to use these tools and do useful things
with them.</p>
<p>Dealing with BLAST can be split up into two steps, both of which can be
done from within Biopython. Firstly, running BLAST for your query
sequence(s), and getting some output. Secondly, parsing the BLAST output
in Python for further analysis.</p>
<p>Your first introduction to running BLAST was probably via the NCBI
web-service. In fact, there are lots of ways you can run BLAST, which
can be categorised several ways. The most important distinction is
running BLAST locally (on your own machine), and running BLAST remotely
(on another machine, typically the NCBI servers). We’re going to start
this chapter by invoking the NCBI online BLAST service from within a
Python script.</p>
<p><em>NOTE</em>: The following Chapter <a class="reference external" href="#chapter:searchio">8</a> describes
<tt class="docutils literal"><span class="pre">Bio.SearchIO</span></tt>, an <em>experimental</em> module in Biopython. We intend this
to ultimately replace the older <tt class="docutils literal"><span class="pre">Bio.Blast</span></tt> module, as it provides a
more general framework handling other related sequence searching tools
as well. However, until that is declared stable, for production code
please continue to use the <tt class="docutils literal"><span class="pre">Bio.Blast</span></tt> module for dealing with NCBI
BLAST.</p>
<div class="section" id="running-blast-over-the-internet">
<h2>7.1  Running BLAST over the Internet<a class="headerlink" href="#running-blast-over-the-internet" title="Permalink to this headline">¶</a></h2>
<p>We use the function <tt class="docutils literal"><span class="pre">qblast()</span></tt> in the <tt class="docutils literal"><span class="pre">Bio.Blast.NCBIWWW</span></tt> module
call the online version of BLAST. This has three non-optional arguments:</p>
<ul class="simple">
<li>The first argument is the blast program to use for the search, as a
lower case string. The options and descriptions of the programs are
available at
<tt class="docutils literal"><span class="pre">`http://www.ncbi.nlm.nih.gov/BLAST/blast_program.shtml</span></tt> &lt;<a class="reference external" href="http://www.ncbi.nlm.nih.gov/BLAST/blast_program.shtml">http://www.ncbi.nlm.nih.gov/BLAST/blast_program.shtml</a>&gt;`__.
Currently <tt class="docutils literal"><span class="pre">qblast</span></tt> only works with blastn, blastp, blastx, tblast
and tblastx.</li>
<li>The second argument specifies the databases to search against. Again,
the options for this are available on the NCBI web pages at
<tt class="docutils literal"><span class="pre">`http://www.ncbi.nlm.nih.gov/BLAST/blast_databases.shtml</span></tt> &lt;<a class="reference external" href="http://www.ncbi.nlm.nih.gov/BLAST/blast_databases.shtml">http://www.ncbi.nlm.nih.gov/BLAST/blast_databases.shtml</a>&gt;`__.</li>
<li>The third argument is a string containing your query sequence. This
can either be the sequence itself, the sequence in fasta format, or
an identifier like a GI number.</li>
</ul>
<p>The <tt class="docutils literal"><span class="pre">qblast</span></tt> function also take a number of other option arguments
which are basically analogous to the different parameters you can set on
the BLAST web page. We’ll just highlight a few of them here:</p>
<ul class="simple">
<li>The <tt class="docutils literal"><span class="pre">qblast</span></tt> function can return the BLAST results in various
formats, which you can choose with the optional <tt class="docutils literal"><span class="pre">format_type</span></tt>
keyword: <tt class="docutils literal"><span class="pre">&quot;HTML&quot;</span></tt>, <tt class="docutils literal"><span class="pre">&quot;Text&quot;</span></tt>, <tt class="docutils literal"><span class="pre">&quot;ASN.1&quot;</span></tt>, or <tt class="docutils literal"><span class="pre">&quot;XML&quot;</span></tt>. The
default is <tt class="docutils literal"><span class="pre">&quot;XML&quot;</span></tt>, as that is the format expected by the parser,
described in section <a class="reference external" href="#sec:parsing-blast">7.3</a> below.</li>
<li>The argument <tt class="docutils literal"><span class="pre">expect</span></tt> sets the expectation or e-value threshold.</li>
</ul>
<p>For more about the optional BLAST arguments, we refer you to the NCBI’s
own documentation, or that built into Biopython:</p>
<p>Note that the default settings on the NCBI BLAST website are not quite
the same as the defaults on QBLAST. If you get different results, you’ll
need to check the parameters (e.g. the expectation value threshold and
the gap values).</p>
<p>For example, if you have a nucleotide sequence you want to search
against the nucleotide database (nt) using BLASTN, and you know the GI
number of your query sequence, you can use:</p>
<p>Alternatively, if we have our query sequence already in a FASTA
formatted file, we just need to open the file and read in this record as
a string, and use that as the query argument:</p>
<p>We could also have read in the FASTA file as a <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> and then
supplied just the sequence itself:</p>
<p>Supplying just the sequence means that BLAST will assign an identifier
for your sequence automatically. You might prefer to use the
<tt class="docutils literal"><span class="pre">SeqRecord</span></tt> object’s format method to make a fasta string (which will
include the existing identifier):</p>
<p>This approach makes more sense if you have your sequence(s) in a
non-FASTA file format which you can extract using <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> (see
Chapter <a class="reference external" href="#chapter:Bio.SeqIO">5</a>).</p>
<p>Whatever arguments you give the <tt class="docutils literal"><span class="pre">qblast()</span></tt> function, you should get
back your results in a handle object (by default in XML format). The
next step would be to parse the XML output into Python objects
representing the search results (Section <a class="reference external" href="#sec:parsing-blast">7.3</a>),
but you might want to save a local copy of the output file first. I find
this especially useful when debugging my code that extracts info from
the BLAST results (because re-running the online search is slow and
wastes the NCBI computer time).</p>
<p>We need to be a bit careful since we can use <tt class="docutils literal"><span class="pre">result_handle.read()</span></tt> to
read the BLAST output only once – calling <tt class="docutils literal"><span class="pre">result_handle.read()</span></tt> again
returns an empty string.</p>
<p>After doing this, the results are in the file <tt class="docutils literal"><span class="pre">my_blast.xml</span></tt> and the
original handle has had all its data extracted (so we closed it).
However, the <tt class="docutils literal"><span class="pre">parse</span></tt> function of the BLAST parser (described
in <a class="reference external" href="#sec:parsing-blast">7.3</a>) takes a file-handle-like object, so
we can just open the saved file for input:</p>
<p>Now that we’ve got the BLAST results back into a handle again, we are
ready to do something with them, so this leads us right into the parsing
section (see Section <a class="reference external" href="#sec:parsing-blast">7.3</a> below). You may want
to jump ahead to that now ….</p>
</div>
<div class="section" id="running-blast-locally">
<h2>7.2  Running BLAST locally<a class="headerlink" href="#running-blast-locally" title="Permalink to this headline">¶</a></h2>
<div class="section" id="introduction">
<h3>7.2.1  Introduction<a class="headerlink" href="#introduction" title="Permalink to this headline">¶</a></h3>
<p>Running BLAST locally (as opposed to over the internet, see
Section <a class="reference external" href="#sec:running-www-blast">7.1</a>) has at least major two
advantages:</p>
<ul class="simple">
<li>Local BLAST may be faster than BLAST over the internet;</li>
<li>Local BLAST allows you to make your own database to search for
sequences against.</li>
</ul>
<p>Dealing with proprietary or unpublished sequence data can be another
reason to run BLAST locally. You may not be allowed to redistribute the
sequences, so submitting them to the NCBI as a BLAST query would not be
an option.</p>
<p>Unfortunately, there are some major drawbacks too – installing all the
bits and getting it setup right takes some effort:</p>
<ul class="simple">
<li>Local BLAST requires command line tools to be installed.</li>
<li>Local BLAST requires (large) BLAST databases to be setup (and
potentially kept up to date).</li>
</ul>
<p>To further confuse matters there are at least four different standalone
BLAST packages, and there are also other tools which can produce
imitation BLAST output files, such as BLAT.</p>
</div>
<div class="section" id="standalone-ncbi-legacy-blast">
<h3>7.2.2  Standalone NCBI “legacy” BLAST<a class="headerlink" href="#standalone-ncbi-legacy-blast" title="Permalink to this headline">¶</a></h3>
<p><a class="reference external" href="http://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&amp;PAGE_TYPE=BlastDocs&amp;DOC_TYPE=Download">NCBI “legacy”
BLAST</a>
included command line tools <tt class="docutils literal"><span class="pre">blastall</span></tt>, <tt class="docutils literal"><span class="pre">blastpgp</span></tt> and <tt class="docutils literal"><span class="pre">rpsblast</span></tt>.
This was the most widely used standalone BLAST tool up until its
replacement BLAST+ was released by the NCBI.</p>
<p>The <tt class="docutils literal"><span class="pre">Bio.Blast.Applications</span></tt> module has wrappers for the “legacy” NCBI
BLAST tools like <tt class="docutils literal"><span class="pre">blastall</span></tt>, <tt class="docutils literal"><span class="pre">blastpgp</span></tt> and <tt class="docutils literal"><span class="pre">rpsblast</span></tt>, and there
are also helper functions in <tt class="docutils literal"><span class="pre">Bio.Blast.NCBIStandalone</span></tt>. These are now
considered obsolete, and will be deprecated and eventually removed from
Biopython as people move over to the replacement BLAST+ suite.</p>
<p>To try and avoid confusion, we will not cover calling these old tools
from Biopython in this tutorial. Have a look at the older edition of
this tutorial included with Biopython 1.52 if you are curious (look at
the Tutorial PDF or HTML file in the Doc directory within
<tt class="docutils literal"><span class="pre">biopython-1.52.tar.gz</span></tt> or <tt class="docutils literal"><span class="pre">biopython-1.52.zip</span></tt>).</p>
</div>
<div class="section" id="standalone-ncbi-blast">
<h3>7.2.3  Standalone NCBI BLAST+<a class="headerlink" href="#standalone-ncbi-blast" title="Permalink to this headline">¶</a></h3>
<p><a class="reference external" href="http://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&amp;PAGE_TYPE=BlastDocs&amp;DOC_TYPE=Download">NCBI “new”
BLAST+</a>
was released in 2009. This replaces the old NCBI “legacy” BLAST package.
The <tt class="docutils literal"><span class="pre">Bio.Blast.Applications</span></tt> module has wrappers for these “new” tools
like <tt class="docutils literal"><span class="pre">blastn</span></tt>, <tt class="docutils literal"><span class="pre">blastp</span></tt>, <tt class="docutils literal"><span class="pre">blastx</span></tt>, <tt class="docutils literal"><span class="pre">tblastn</span></tt>, <tt class="docutils literal"><span class="pre">tblastx</span></tt> (which
all used to be handled by <tt class="docutils literal"><span class="pre">blastall</span></tt>), <tt class="docutils literal"><span class="pre">psiblast</span></tt> (replacing
<tt class="docutils literal"><span class="pre">blastpgp</span></tt>) and <tt class="docutils literal"><span class="pre">rpsblast</span></tt> and <tt class="docutils literal"><span class="pre">rpstblastn</span></tt> (which replace the old
<tt class="docutils literal"><span class="pre">rpsblast</span></tt>). We don’t include a wrapper for the <tt class="docutils literal"><span class="pre">makeblastdb</span></tt> used
in BLAST+ to build a local BLAST database from FASTA file, nor the
equivalent tool <tt class="docutils literal"><span class="pre">formatdb</span></tt> in “legacy” BLAST.</p>
<p>This section will show briefly how to use these tools from within
Python. If you have already read or tried the alignment tool examples in
Section <a class="reference external" href="#sec:alignment-tools">6.4</a> this should all seem quite
straightforward. First, we construct a command line string (as you would
type in at the command line prompt if running standalone BLAST by hand).
Then we can execute this command from within Python.</p>
<p>For example, taking a FASTA file of gene nucleotide sequences, you might
want to run a BLASTX (translation) search against the non-redundant (NR)
protein database. Assuming you (or your systems administrator) has
downloaded and installed the NR database, you might run:</p>
<p>This should run BLASTX against the NR database, using an expectation
cut-off value of 0.001 and produce XML output to the specified file
(which we can then parse). On my computer this takes about six minutes -
a good reason to save the output to a file so you and repeat any
analysis as needed.</p>
<p>From within Biopython we can use the NCBI BLASTX wrapper from the
<tt class="docutils literal"><span class="pre">Bio.Blast.Applications</span></tt> module to build the command line string, and
run it:</p>
<p>In this example there shouldn’t be any output from BLASTX to the
terminal, so stdout and stderr should be empty. You may want to check
the output file <tt class="docutils literal"><span class="pre">opuntia.xml</span></tt> has been created.</p>
<p>As you may recall from earlier examples in the tutorial, the
<tt class="docutils literal"><span class="pre">opuntia.fasta</span></tt> contains seven sequences, so the BLAST XML output
should contain multiple results. Therefore use
<tt class="docutils literal"><span class="pre">Bio.Blast.NCBIXML.parse()</span></tt> to parse it as described below in
Section <a class="reference external" href="#sec:parsing-blast">7.3</a>.</p>
</div>
<div class="section" id="wu-blast-and-ab-blast">
<h3>7.2.4  WU-BLAST and AB-BLAST<a class="headerlink" href="#wu-blast-and-ab-blast" title="Permalink to this headline">¶</a></h3>
<p>You may also come across <a class="reference external" href="http://blast.wustl.edu/">Washington University
BLAST</a> (WU-BLAST), and its successor,
<a class="reference external" href="http://blast.advbiocomp.com">Advanced Biocomputing BLAST</a> (AB-BLAST,
released in 2009, not free/open source). These packages include the
command line tools <tt class="docutils literal"><span class="pre">wu-blastall</span></tt> and <tt class="docutils literal"><span class="pre">ab-blastall</span></tt>.</p>
<p>Biopython does not currently provide wrappers for calling these tools,
but should be able to parse any NCBI compatible output from them.</p>
</div>
</div>
<div class="section" id="parsing-blast-output">
<h2>7.3  Parsing BLAST output<a class="headerlink" href="#parsing-blast-output" title="Permalink to this headline">¶</a></h2>
<p>As mentioned above, BLAST can generate output in various formats, such
as XML, HTML, and plain text. Originally, Biopython had parsers for
BLAST plain text and HTML output, as these were the only output formats
offered at the time. Unfortunately, the BLAST output in these formats
kept changing, each time breaking the Biopython parsers. Our HTML BLAST
parser has been removed, but the plain text BLAST parser is still
available (see Section <a class="reference external" href="#sec:parsing-blast-deprecated">7.5</a>). Use
it at your own risk, it may or may not work, depending on which BLAST
version you’re using.</p>
<p>As keeping up with changes in BLAST became a hopeless endeavor,
especially with users running different BLAST versions, we now recommend
to parse the output in XML format, which can be generated by recent
versions of BLAST. Not only is the XML output more stable than the plain
text and HTML output, it is also much easier to parse automatically,
making Biopython a whole lot more stable.</p>
<p>You can get BLAST output in XML format in various ways. For the parser,
it doesn’t matter how the output was generated, as long as it is in the
XML format.</p>
<ul class="simple">
<li>You can use Biopython to run BLAST over the internet, as described in
section <a class="reference external" href="#sec:running-www-blast">7.1</a>.</li>
<li>You can use Biopython to run BLAST locally, as described in
section <a class="reference external" href="#sec:running-local-blast">7.2</a>.</li>
<li>You can do the BLAST search yourself on the NCBI site through your
web browser, and then save the results. You need to choose XML as the
format in which to receive the results, and save the final BLAST page
you get (you know, the one with all of the interesting results!) to a
file.</li>
<li>You can also run BLAST locally without using Biopython, and save the
output in a file. Again, you need to choose XML as the format in
which to receive the results.</li>
</ul>
<p>The important point is that you do not have to use Biopython scripts to
fetch the data in order to be able to parse it. Doing things in one of
these ways, you then need to get a handle to the results. In Python, a
handle is just a nice general way of describing input to any info source
so that the info can be retrieved using <tt class="docutils literal"><span class="pre">read()</span></tt> and <tt class="docutils literal"><span class="pre">readline()</span></tt>
functions (see Section sec:appendix-handles).</p>
<p>If you followed the code above for interacting with BLAST through a
script, then you already have <tt class="docutils literal"><span class="pre">result_handle</span></tt>, the handle to the BLAST
results. For example, using a GI number to do an online search:</p>
<p>If instead you ran BLAST some other way, and have the BLAST output (in
XML format) in the file <tt class="docutils literal"><span class="pre">my_blast.xml</span></tt>, all you need to do is to open
the file for reading:</p>
<p>Now that we’ve got a handle, we are ready to parse the output. The code
to parse it is really quite small. If you expect a single BLAST result
(i.e. you used a single query):</p>
<p>or, if you have lots of results (i.e. multiple query sequences):</p>
<p>Just like <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> and <tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt> (see
Chapters <a class="reference external" href="#chapter:Bio.SeqIO">5</a>
and <a class="reference external" href="#chapter:Bio.AlignIO">6</a>), we have a pair of input functions,
<tt class="docutils literal"><span class="pre">read</span></tt> and <tt class="docutils literal"><span class="pre">parse</span></tt>, where <tt class="docutils literal"><span class="pre">read</span></tt> is for when you have exactly one
object, and <tt class="docutils literal"><span class="pre">parse</span></tt> is an iterator for when you can have lots of
objects – but instead of getting <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> or
<tt class="docutils literal"><span class="pre">MultipleSeqAlignment</span></tt> objects, we get BLAST record objects.</p>
<p>To be able to handle the situation where the BLAST file may be huge,
containing thousands of results, <tt class="docutils literal"><span class="pre">NCBIXML.parse()</span></tt> returns an
iterator. In plain English, an iterator allows you to step through the
BLAST output, retrieving BLAST records one by one for each BLAST search
result:</p>
<p>Or, you can use a <tt class="docutils literal"><span class="pre">for</span></tt>-loop:</p>
<p>Note though that you can step through the BLAST records only once.
Usually, from each BLAST record you would save the information that you
are interested in. If you want to save all returned BLAST records, you
can convert the iterator into a list:</p>
<p>Now you can access each BLAST record in the list with an index as usual.
If your BLAST file is huge though, you may run into memory problems
trying to save them all in a list.</p>
<p>Usually, you’ll be running one BLAST search at a time. Then, all you
need to do is to pick up the first (and only) BLAST record in
<tt class="docutils literal"><span class="pre">blast_records</span></tt>:</p>
<p>or more elegantly:</p>
<p>I guess by now you’re wondering what is in a BLAST record.</p>
</div>
<div class="section" id="the-blast-record-class">
<h2>7.4  The BLAST record class<a class="headerlink" href="#the-blast-record-class" title="Permalink to this headline">¶</a></h2>
<p>A BLAST Record contains everything you might ever want to extract from
the BLAST output. Right now we’ll just show an example of how to get
some info out of the BLAST report, but if you want something in
particular that is not described here, look at the info on the record
class in detail, and take a gander into the code or automatically
generated documentation – the docstrings have lots of good info about
what is stored in each piece of information.</p>
<p>To continue with our example, let’s just print out some summary info
about all hits in our blast report greater than a particular threshold.
The following code does this:</p>
<p>This will print out summary reports like the following:</p>
<p>Basically, you can do anything you want to with the info in the BLAST
report once you have parsed it. This will, of course, depend on what you
want to use it for, but hopefully this helps you get started on doing
what you need to do!</p>
<p>An important consideration for extracting information from a BLAST
report is the type of objects that the information is stored in. In
Biopython, the parsers return <tt class="docutils literal"><span class="pre">Record</span></tt> objects, either <tt class="docutils literal"><span class="pre">Blast</span></tt> or
<tt class="docutils literal"><span class="pre">PSIBlast</span></tt> depending on what you are parsing. These objects are
defined in <tt class="docutils literal"><span class="pre">Bio.Blast.Record</span></tt> and are quite complete.</p>
<p>Here are my attempts at UML class diagrams for the <tt class="docutils literal"><span class="pre">Blast</span></tt> and
<tt class="docutils literal"><span class="pre">PSIBlast</span></tt> record classes. If you are good at UML and see
mistakes/improvements that can be made, please let me know. The Blast
class diagram is shown in Figure <a class="reference external" href="#fig:blastrecord">7.4</a>.</p>
<p><a href="#id1"><span class="problematic" id="id2">|image1|</span></a></p>
<p>The PSIBlast record object is similar, but has support for the rounds
that are used in the iteration steps of PSIBlast. The class diagram for
PSIBlast is shown in Figure <a class="reference external" href="#fig:psiblastrecord">7.4</a>.</p>
<p><a href="#id3"><span class="problematic" id="id4">|image2|</span></a></p>
</div>
<div class="section" id="deprecated-blast-parsers">
<h2>7.5  Deprecated BLAST parsers<a class="headerlink" href="#deprecated-blast-parsers" title="Permalink to this headline">¶</a></h2>
<p>Older versions of Biopython had parsers for BLAST output in plain text
or HTML format. Over the years, we discovered that it is very hard to
maintain these parsers in working order. Basically, any small change to
the BLAST output in newly released BLAST versions tends to cause the
plain text and HTML parsers to break. We therefore recommend parsing
BLAST output in XML format, as described in
section <a class="reference external" href="#sec:parsing-blast">7.3</a>.</p>
<p>Depending on which BLAST versions or programs you’re using, our plain
text BLAST parser may or may not work. Use it at your own risk!</p>
<div class="section" id="parsing-plain-text-blast-output">
<h3>7.5.1  Parsing plain-text BLAST output<a class="headerlink" href="#parsing-plain-text-blast-output" title="Permalink to this headline">¶</a></h3>
<p>The plain text BLAST parser is located in <tt class="docutils literal"><span class="pre">Bio.Blast.NCBIStandalone</span></tt>.</p>
<p>As with the XML parser, we need to have a handle object that we can pass
to the parser. The handle must implement the <tt class="docutils literal"><span class="pre">readline()</span></tt> method and
do this properly. The common ways to get such a handle are to either use
the provided <tt class="docutils literal"><span class="pre">blastall</span></tt> or <tt class="docutils literal"><span class="pre">blastpgp</span></tt> functions to run the local
blast, or to run a local blast via the command line, and then do
something like the following:</p>
<p>Well, now that we’ve got a handle (which we’ll call <tt class="docutils literal"><span class="pre">result_handle</span></tt>),
we are ready to parse it. This can be done with the following code:</p>
<p>This will parse the BLAST report into a Blast Record class (either a
Blast or a PSIBlast record, depending on what you are parsing) so that
you can extract the information from it. In our case, let’s just use
print out a quick summary of all of the alignments greater than some
threshold value.</p>
<p>If you also read the section <a class="reference external" href="#sec:parsing-blast">7.3</a> on parsing
BLAST XML output, you’ll notice that the above code is identical to what
is found in that section. Once you parse something into a record class
you can deal with it independent of the format of the original BLAST
info you were parsing. Pretty snazzy!</p>
<p>Sure, parsing one record is great, but I’ve got a BLAST file with tons
of records – how can I parse them all? Well, fear not, the answer lies
in the very next section.</p>
</div>
<div class="section" id="parsing-a-plain-text-blast-file-full-of-blast-runs">
<h3>7.5.2  Parsing a plain-text BLAST file full of BLAST runs<a class="headerlink" href="#parsing-a-plain-text-blast-file-full-of-blast-runs" title="Permalink to this headline">¶</a></h3>
<p>We can do this using the blast iterator. To set up an iterator, we first
set up a parser, to parse our blast reports in Blast Record objects:</p>
<p>Then we will assume we have a handle to a bunch of blast records, which
we’ll call <tt class="docutils literal"><span class="pre">result_handle</span></tt>. Getting a handle is described in full
detail above in the blast parsing sections.</p>
<p>Now that we’ve got a parser and a handle, we are ready to set up the
iterator with the following command:</p>
<p>The second option, the parser, is optional. If we don’t supply a parser,
then the iterator will just return the raw BLAST reports one at a time.</p>
<p>Now that we’ve got an iterator, we start retrieving blast records
(generated by our parser) using <tt class="docutils literal"><span class="pre">next()</span></tt>:</p>
<p>Each call to next will return a new record that we can deal with. Now we
can iterate through this records and generate our old favorite, a nice
little blast report:</p>
<p>The iterator allows you to deal with huge blast records without any
memory problems, since things are read in one at a time. I have parsed
tremendously huge files without any problems using this.</p>
</div>
<div class="section" id="finding-a-bad-record-somewhere-in-a-huge-plain-text-blast-file">
<h3>7.5.3  Finding a bad record somewhere in a huge plain-text BLAST file<a class="headerlink" href="#finding-a-bad-record-somewhere-in-a-huge-plain-text-blast-file" title="Permalink to this headline">¶</a></h3>
<p>One really ugly problem that happens to me is that I’ll be parsing a
huge blast file for a while, and the parser will bomb out with a
ValueError. This is a serious problem, since you can’t tell if the
ValueError is due to a parser problem, or a problem with the BLAST. To
make it even worse, you have no idea where the parse failed, so you
can’t just ignore the error, since this could be ignoring an important
data point.</p>
<p>We used to have to make a little script to get around this problem, but
the <tt class="docutils literal"><span class="pre">Bio.Blast</span></tt> module now includes a <tt class="docutils literal"><span class="pre">BlastErrorParser</span></tt> which
really helps make this easier. The <tt class="docutils literal"><span class="pre">BlastErrorParser</span></tt> works very
similar to the regular <tt class="docutils literal"><span class="pre">BlastParser</span></tt>, but it adds an extra layer of
work by catching ValueErrors that are generated by the parser, and
attempting to diagnose the errors.</p>
<p>Let’s take a look at using this parser – first we define the file we are
going to parse and the file to write the problem reports to:</p>
<p>Now we want to get a <tt class="docutils literal"><span class="pre">BlastErrorParser</span></tt>:</p>
<p>Notice that the parser take an optional argument of a handle. If a
handle is passed, then the parser will write any blast records which
generate a ValueError to this handle. Otherwise, these records will not
be recorded.</p>
<p>Now we can use the <tt class="docutils literal"><span class="pre">BlastErrorParser</span></tt> just like a regular blast
parser. Specifically, we might want to make an iterator that goes
through our blast records one at a time and parses them with the error
parser:</p>
<p>We can read these records one a time, but now we can catch and deal with
errors that are due to problems with Blast (and not with the parser
itself):</p>
<p>The <tt class="docutils literal"><span class="pre">.next()</span></tt> method is normally called indirectly via a <tt class="docutils literal"><span class="pre">for</span></tt>-loop.
Right now the <tt class="docutils literal"><span class="pre">BlastErrorParser</span></tt> can generate the following errors:</p>
<ul class="simple">
<li><tt class="docutils literal"><span class="pre">ValueError</span></tt> – This is the same error generated by the regular
BlastParser, and is due to the parser not being able to parse a
specific file. This is normally either due to a bug in the parser, or
some kind of discrepancy between the version of BLAST you are using
and the versions the parser is able to handle.</li>
<li><tt class="docutils literal"><span class="pre">LowQualityBlastError</span></tt> – When BLASTing a sequence that is of really
bad quality (for example, a short sequence that is basically a
stretch of one nucleotide), it seems that Blast ends up masking out
the entire sequence and ending up with nothing to parse. In this case
it will produce a truncated report that causes the parser to generate
a ValueError. <tt class="docutils literal"><span class="pre">LowQualityBlastError</span></tt> is reported in these cases.
This error returns an info item with the following information:<ul>
<li><tt class="docutils literal"><span class="pre">item[0]</span></tt> – The error message</li>
<li><tt class="docutils literal"><span class="pre">item[1]</span></tt> – The id of the input record that caused the error.
This is really useful if you want to record all of the records
that are causing problems.</li>
</ul>
</li>
</ul>
<p>As mentioned, with each error generated, the BlastErrorParser will write
the offending record to the specified <tt class="docutils literal"><span class="pre">error_handle</span></tt>. You can then go
ahead and look and these and deal with them as you see fit. Either you
will be able to debug the parser with a single blast report, or will
find out problems in your blast runs. Either way, it will definitely be
a useful experience!</p>
<p>Hopefully the <tt class="docutils literal"><span class="pre">BlastErrorParser</span></tt> will make it much easier to debug and
deal with large Blast files.</p>
</div>
</div>
<div class="section" id="dealing-with-psi-blast">
<h2>7.6  Dealing with PSI-BLAST<a class="headerlink" href="#dealing-with-psi-blast" title="Permalink to this headline">¶</a></h2>
<p>You can run the standalone version of PSI-BLAST (the legacy NCBI command
line tool <tt class="docutils literal"><span class="pre">blastpgp</span></tt>, or its replacement <tt class="docutils literal"><span class="pre">psiblast</span></tt>) using the
wrappers in <tt class="docutils literal"><span class="pre">Bio.Blast.Applications</span></tt> module.</p>
<p>At the time of writing, the NCBI do not appear to support tools running
a PSI-BLAST search via the internet.</p>
<p>Note that the <tt class="docutils literal"><span class="pre">Bio.Blast.NCBIXML</span></tt> parser can read the XML output from
current versions of PSI-BLAST, but information like which sequences in
each iteration is new or reused isn’t present in the XML file. If you
care about this information you may have more joy with the plain text
output and the <tt class="docutils literal"><span class="pre">PSIBlastParser</span></tt> in <tt class="docutils literal"><span class="pre">Bio.Blast.NCBIStandalone</span></tt>.</p>
</div>
<div class="section" id="dealing-with-rps-blast">
<h2>7.7  Dealing with RPS-BLAST<a class="headerlink" href="#dealing-with-rps-blast" title="Permalink to this headline">¶</a></h2>
<p>You can run the standalone version of RPS-BLAST (either the legacy NCBI
command line tool <tt class="docutils literal"><span class="pre">rpsblast</span></tt>, or its replacement with the same name)
using the wrappers in <tt class="docutils literal"><span class="pre">Bio.Blast.Applications</span></tt> module.</p>
<p>At the time of writing, the NCBI do not appear to support tools running
an RPS-BLAST search via the internet.</p>
<p>You can use the <tt class="docutils literal"><span class="pre">Bio.Blast.NCBIXML</span></tt> parser to read the XML output from
current versions of RPS-BLAST.</p>
</div>
</div>


          </div>
        </div>
      </div>
      <div class="sphinxsidebar">
        <div class="sphinxsidebarwrapper">
  <h3><a href="index.html">Table Of Contents</a></h3>
  <ul>
<li><a class="reference internal" href="#">Chapter 7  BLAST</a><ul>
<li><a class="reference internal" href="#running-blast-over-the-internet">7.1  Running BLAST over the Internet</a></li>
<li><a class="reference internal" href="#running-blast-locally">7.2  Running BLAST locally</a><ul>
<li><a class="reference internal" href="#introduction">7.2.1  Introduction</a></li>
<li><a class="reference internal" href="#standalone-ncbi-legacy-blast">7.2.2  Standalone NCBI “legacy” BLAST</a></li>
<li><a class="reference internal" href="#standalone-ncbi-blast">7.2.3  Standalone NCBI BLAST+</a></li>
<li><a class="reference internal" href="#wu-blast-and-ab-blast">7.2.4  WU-BLAST and AB-BLAST</a></li>
</ul>
</li>
<li><a class="reference internal" href="#parsing-blast-output">7.3  Parsing BLAST output</a></li>
<li><a class="reference internal" href="#the-blast-record-class">7.4  The BLAST record class</a></li>
<li><a class="reference internal" href="#deprecated-blast-parsers">7.5  Deprecated BLAST parsers</a><ul>
<li><a class="reference internal" href="#parsing-plain-text-blast-output">7.5.1  Parsing plain-text BLAST output</a></li>
<li><a class="reference internal" href="#parsing-a-plain-text-blast-file-full-of-blast-runs">7.5.2  Parsing a plain-text BLAST file full of BLAST runs</a></li>
<li><a class="reference internal" href="#finding-a-bad-record-somewhere-in-a-huge-plain-text-blast-file">7.5.3  Finding a bad record somewhere in a huge plain-text BLAST file</a></li>
</ul>
</li>
<li><a class="reference internal" href="#dealing-with-psi-blast">7.6  Dealing with PSI-BLAST</a></li>
<li><a class="reference internal" href="#dealing-with-rps-blast">7.7  Dealing with RPS-BLAST</a></li>
</ul>
</li>
</ul>

  <h4>Previous topic</h4>
  <p class="topless"><a href="chr6.html"
                        title="previous chapter">Chapter 6  Multiple Sequence Alignment objects</a></p>
  <h4>Next topic</h4>
  <p class="topless"><a href="chr8.html"
                        title="next chapter">第8章  BLAST和其他序列搜索工具(<em>实验性质的代码</em>)</a></p>
  <h3>This Page</h3>
  <ul class="this-page-menu">
    <li><a href="_sources/chr7.txt"
           rel="nofollow">Show Source</a></li>
  </ul>
<div id="searchbox" style="display: none">
  <h3>Quick search</h3>
    <form class="search" action="search.html" method="get">
      <input type="text" name="q" />
      <input type="submit" value="Go" />
      <input type="hidden" name="check_keywords" value="yes" />
      <input type="hidden" name="area" value="default" />
    </form>
    <p class="searchtip" style="font-size: 90%">
    Enter search terms or a module, class or function name.
    </p>
</div>
<script type="text/javascript">$('#searchbox').show(0);</script>
        </div>
      </div>
      <div class="clearer"></div>
    </div>
    <div class="related">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="genindex.html" title="General Index"
             >index</a></li>
        <li class="right" >
          <a href="chr8.html" title="第8章 BLAST和其他序列搜索工具(实验性质的代码)"
             >next</a> |</li>
        <li class="right" >
          <a href="chr6.html" title="Chapter 6 Multiple Sequence Alignment objects"
             >previous</a> |</li>
        <li><a href="index.html">Biopython_en 1.0 documentation</a> &raquo;</li> 
      </ul>
    </div>
    <div class="footer">
        &copy; Copyright 2013, Biopython.
      Created using <a href="http://sphinx-doc.org/">Sphinx</a> 1.2b1.
    </div>
  </body>
</html>